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Calculations of formation energies and charge transition levels of defects routinely rely on density 
functional theory (DFT) for describing the electronic structure. Since bulk band gaps of semi- 
conductors and insulators are not well described in semilocal approximations to DFT, band-gap 
correction schemes or advanced theoretical models which properly describe band gaps need to be 
employed. However, it has become apparent that different methods that reproduce the experimental 
band gap can yield substantially different results regarding charge transition levels of point defects. 
We investigate this problem in the case of the (-1-2/0) charge transition level of the O vacancy in 
ZnO, which has attracted considerable attention as a benchmark case. For this purpose, we first 
perform calculations based on non-screened hybrid density functionals, and then compare our results 
with those of other methods. While our results agree very well with those obtained with screened 
hybrid functionals, they are strikingly different compared to those obtained with other band-gap- 
corrected schemes. Nevertheless, we show that all the different methods agree well with each other 
and with our calculations when a suitable alignment procedure is adopted. The proposed procedure 
consists in aligning the electron band structure through an external potential, such as the vacuum 
level. When the electron densities are well reproduced, this procedure is equivalent to an alignment 
through the average electrostatic potential in a calculation subject to periodic boundary conditions. 
We stress that, in order to give accurate defect levels, a theoretical scheme is required to yield not 
only band gaps in agreement with experiment, but also band edges correctly positioned with respect 
to such a reference potential. 

PACS numbers: 71.15.Nc, 71.55.-i, 71.55.Gs 



I. INTRODUCTION 

Point defects can affect the properties of solids in a 
dramatic way,^ They determine, for example, the con- 
ductivity of semiconductors, the color of natural crys- 
tals, and the mechanical properties of materials. Equally 
important, defects influence or govern the performance 
and the long-term stability of a wide range of semicon- 
ductor devices, such as metal-oxide-semiconductor field- 
effect transistors, photovoltaic cells, solid fuel cells, to 
name a few. The theoretical characterization of de- 
fects, especially in wide band-gap materials, has be- 
come increasingly important in the attempt to under- 
stand and control the performance of these devices 
In the last decades, density functional theory (DFT) 
has grown into the standard theoretical model to de- 
scribe the electronic and atomic structure of solids. The 
common approximations to DFT, viz. the local density 
approximation (LDA) and the generalized gradient ap- 
proximation (GGA), systematically underestimate band 
gaps of semiconductors and insulators. Since the band 
gap is the relevant energy scale in the study of defects, 
this so-called "band-gap problem" of LDA and GGA 
severely affects the predictive power of these approxi- 
mations when applied to defect levels. Recently there 
have been lots of efforts to assess the importance of band 
gap corrections^"'' and to use theoretical models giving 
a much more appropriate description of the bulk band 
structure. The choice of methods is large and includes the 



LDA+U method^'"— approximate self- interaction cor- 
rection schemeSfi^ hybrid density functionals^^"— the 
use of modified pseudopotentials,-- empirical schemes^^ 
and more advanced theoretical tools, such as the many- 
body perturbation theory within the GW and higher 
approximations ,21^ 

It appears evident to assume that a good theoretical 
model must at least satisfy two conditions, namely (i) 
give an accurate electron density of the defect system and 
(ii) yield a good band gap of the host material. While 
these two requirements form a necessary prerequisite to 
obtain reliable results concerning defect formation ener- 
gies and associated charge transition levelSj— "4 it has re- 
cently become apparent that it is by no means sufficient. 
This is best exemplified in the case of defect energy levels 
in ZnOj^i^"— lii'Sii^"-^^ This is a particularly severe case, 
because the LDA and the GGA yield a bulk band-gap 
of 0.6-0.8 eV, severely underestimating the experimen- 
tal value of 3.44 eV. For the case of the (+2/0) charge 
transition level of the oxygen vacancy (Vo) theoretical 
models yield levels either as low as 0.6 eV above the va- 
lence band maximum (VBM) or as high as 2.4 eV above 
VBM. These results differ significantly despite the fact 
that in all these theoretical models the "band-gap prob- 
lem" was accounted for. In addition, other critical issues, 
such as finite-size effects associated to the supercell treat- 
ment, were presumably under control in these studies. 
Furthermore, the first condition concerning the accuracy 
of the electron density was clearly also fulfilled since the 



2 



involved electronic state corresponds to the fully sym- 
metric Oi state which is already correctly described via 
a semilocal functional. The second condition concerning 
the band gap was fulfilled by construction. 

Recently, Lany and Zunger provided a very detailed 
overview of the way various theoretical and computa- 
tional approximations affect the determination of defect 
formation energies and charge transition levels.-^' They 
concluded that, in addition to the two requirements dis- 
cussed above, a reliable theoretical model should cor- 
rectly describe the relative positions of all relevant elec- 
tronic states. For ZnO, this condition mainly concerns 
the position of the Zn 3d states with respect to the 
conduction and valence band edges. The importance 
of this requirement becomes evident when considering 
shallow defects, the wavefunctions of which can be al- 
ways thought as arising from a linear combination of bulk 
bands. 

In this work, we show that that there is yet another 
crucial requirement that the theoretical model must ful- 
fill. In order to yield an appropriate description of defect 
formation energies and associated charge transition lev- 
els, the positions of the VBM and the conduction band 
minimum (CBM) with respect to a suitably defined ref- 
erence potential should also be accurately described. To 
demonstrate this, we first calculate the (+2/0) charge 
transition level of the Vq in ZnO and compare our result 
with those available in the literature. Our study adds 
to a series of studies in which conflicting re- 

sults were found. However, we show that these seem- 
ingly incompatible findings agree reasonably with each 
other when an alternative alignment scheme is used. We 
provide theoretical arguments to rationalize this finding. 
Similar results are expected for other atomically localized 
defects and for other materials in which the "band-gap 
problem" of semilocal calculations is particularly severe. 
Our investigation thus leads to a deeper understanding of 
the "band-edge problem" in the theoretical study of de- 
fect levels and provides a requirement for the theoretical 
model in addition to the conditions mentioned above. 

This paper is organized as follows. In Sec. |TT1 we sum- 
marize our computational approach for calculating defect 
formation energies and charge transition levels. The ob- 
tained results are discussed and compared to other cal- 
culations in Sec. IIIII An alignment scheme with respect 
to the average electrostatic potential is introduced and 
found to bring all the calculated results in good agree- 
ment with each other. The significance of this alignment 
of bulk band structures is discussed in more detail in Sec. 
IIVI To understand our findings about defect charge tran- 
sition levels, fundamental differences between localized 
and extended states in approximate DFT formulations 
are discussed in Sec. [V] In Sec. |Vll two different theo- 
ries reproducing the experimental band gap but differing 
in the positions of the bulk band edges with respect to 
the vacuum level are taken under consideration to com- 
plete our rationale. We summarize our work and draw 
conclusions in Sec. IVIII 



II. COMPUTATIONAL METHODS 

In the present calculations, the electronic structure was 
treated using two different functionals. First, we em- 
ployed the GGA functional proposed by Perdew, Burke, 
and Ernherhof (PBE).'^^ For comparison with previous 
calculations in the literature, we obtained for bulk ZnO 
a band gap of 0.83 eV, to be compared with the exper- 
imental value of 3.44 eV. To obtain an improved band 
gap, we used a hybrid density functional'^- defined by a 
single parameter a corresponding to the fraction of non- 
local Fock exchange admixed to the GGA exchange: 

£,hybrid^^^Fock^(^_^)^GGA^ (1) 

A hybrid functional with a — 0.25 and with the PBE 
for the GGA part^ is referred to as PBEO, PBEh, or 
PBEIPBE. For ZnO, we obtained a band gap of 2.82 eV 
using this functional. The experimental band gap is re- 
produced with a = 0.32. In the following, we refer to 
this functional as to PBEh-32. While this adjustment of 
a is empirical, it can be justified to a certain extent.—"— 
It can be shown that the optimal value of Oopt, i.e. the 
one which reproduces the experimental band gap, is ap- 
proximately given by Oopt ^ l/^oo- Here, £oo is the elec- 
tronic part of the static dielectric constant. For a large 
number of materials this relationship is approximately 
fulfilled4Si^ The adjustment of a can also be justified 
in some cases by comparison with more accurate GW 
calculations."^^ 

The main quantity that needs to be calculated is the 
formation energy of the oxygen vacancy in a charge state 
q, which is given as)^ 

= E^^, - £:tot,buik + MO + q{ev + £f)- (2) 

Here E^^^. is the total energy of the defect system con- 
taining a single O vacancy in the supercell, i^tot.buik is 
the total energy of the host material without any defect, 
fiQ is the atomic chemical potential of oxygen, and s-p is 
the electron chemical potential. The latter is referred to 
the VBM ey. Except for semiconductors with degener- 
ate doping, e-p varies between zero and the band gap of 
the material E„. 

The atomic chemical potentials /io a-nd /izn are bound 
by the condition that ZnO is in thermal equilibrium with 
the reservoir of O and Zn atoms, i.e. ^Zn + Mo = MZnO- 
Oxygen-rich conditions are defined by the onset of spon- 
taneous formation of O2 molecules, i.e. by /io — ^E^^l. 
Oxygen-poor (Zn-rich) conditions are correspondingly 
defined by the onset of spontaneous formation of bulk 
Zn crystallites, i.e. via ^zn = Mzn,buik- The formation of 
oxygen vacancies in ZnO is hindered in 0-rich conditions, 
and facilitated in 0-poor conditions. The calculation of 
the O chemical potential in 0-poor conditions poses some 
difficulties when hybrid density functionals are used, be- 
cause this involves the calculation of the total energy of 
bulk Zn. In Hartree-Fock theory the description of metals 
leads to divergences, and the same problem is also found 
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with hybrid functionals. To overcome this problem, we 
assume that the cohesive energy of bulk Zn, which is 
well described in the GGA, does not change significantly 
in the hybrid functional calculation.^^ Alternatively, one 
could define the O chemical potential in 0-poor condi- 
tions by assuming that the separation between the 0-rich 
and 0-poor chemical potentials in GGA is preserved in 
the hybrid functional calculation; this condition corre- 
sponds to assuming equal formation energies for ZnO in 
GGA and in the hybrid functional scheme. These two 
ways of determining the O chemical potential in 0-poor 
conditions lead to formation energies differing by about 
0.4 eV. 

Gharge transition levels correspond to the specific 
value of the electron chemical potential for which two 
charge states have equal formation energies. The (-1-2/0) 
charge transition level is thus given by: 



-2/0) 



+2 



■ey. 



(3) 



Charge transition levels do not depend on atomic chem- 
ical potentials. 

The calculations were performed within a plane- 
wave pseudopotcntial formulation. Soft norm-conserving 
pseudopotentials^ were generated at the PBE level and 
used in all subsequent calculations. The plane- wave 
kinetic energy cutoff, determined by the much harder 
O pseudopotcntial, was set to 80 Ry. The calcula- 
tions in the present paper were performed with the code 
CPMd4^"— We explicitly treated the singularity in the 
nonlocal exchange potential. 

We used the experimental lattice parameters for bulk 
ZnO, since these were found to be very close to theoreti- 
cal lattice parameters obtained with hybrid functionals 
We also used experimental lattice constants in our GGA 
calculations, finding results which did not differ in any 
significant way from GGA calculations performed with 
theoretical lattice parameters^ Upon defect formation, 
geometry relaxations were performed with both the GGA 
and the PBEh-32 functionals. The defect structures 
achieved in the two cases were found to be very similar: in 
PBEh-32, for example, PBE-optimized defect structures 
are only 0.08 eV higher in energy than those optimized 
consistently at the PBEh-32 level. Hence, geometry op- 
timization at the PBEh-32 level has no effect on the po- 
sition of the (-1-2/0) charge transition level [Eq. 

For the defect structures we used the supercell ap- 
proach. This gives rise to finite-size effects which need 
to be accounted for. First, as suggested by Van de Walle 
and Neugebauer,^ the total energies of charged defects 
were corrected by qAV, AV being the shift needed to 
align the local potential of the neutral system far from 
the defect to that of a separate unperturbed bulk calcu- 
lation, which was used to determine sy. This term was 
found to be quite small for the supercells employed in 
our calculations. Second, the total energies of charged 
defect states are subject to spurious electrostatic contri- 
butions associated to the periodic boundary conditions 
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FIG. 1; (Color online) Charge transition level £(+2/0) vs 
inverse number of atoms contained in the supercell A^at, (a) 
for the PBE calculation (1x1x1 and 2x2x2 fc-point meshes) 
and (b) for the PBEh-32 calculation (1x1x1 mesh). £ (-1-2/0) 
is referred to the respective VBM. 



and to the compensating background charge in our su- 
percell calculations. To evaluate these effects, we used 
an extrapolation scheme based on supercell calculations 
of increasing size, containing 96, 192, and 384 atoms, as 
shown in Fig. [TJ When using the PBE functional, the 
convergence of formation energies and charge transition 
levels is accelerated when using the 2x2x2 Monkhorst- 
Pack mesh instead of a sampling at the sole F point [Fig. 
[IJa)]. Hence, finite-size corrections are sizeable for the 
PBE calculation and a careful extrapolation of the re- 
sults is needed, as previously shown by Oba et al^ At 
variance, a denser /c-point mesh turned out to be unneces- 
sary for a calculation with the hybrid functional PBEh-32 
[Fig. mb)]. Indeed, in the latter case, the bulk band gap 
is substantially larger and the dispersion of the defect 
state is already negligible for the smallest supercells con- 
sidered. This behavior is in line with observations in a 
previous study on defects in ZnO.''^ A notable difference 
between finite size effects in PBE and PBEh-32 calcula- 
tions suggests that unphysical defect-defect interactions 
mediated by bulk bands could be operative in the former 
caseA For the largest supercell considered here, we ob- 
tain a conservative estimate of 0.20 eV for the residual 
finite-size error by considering the monopole correction 
proposed by Makov and Payne4^ 



III. OXYGEN VACANCY IN ZnO 

For the neutral oxygen vacancy, we obtained, at the 
PBE level, formation energies of 3.17 eV in 0-rich condi- 
tions and of 0.50 eV in 0-poor conditions. In the PBEh- 
32 calculation, the corresponding value is 3.57 eV in O- 
rich conditions. In 0-poor conditions, we found 0.50 and 
0.90 eV depending on whether the cohesive energy of Zn 
or the formation energy of ZnO is taken from the GGA, 
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FIG. 2: (Color online) Formation energy of oxygen vacancy 
in ZnO vs electron chemical potential, as obtained with the 
PBEh-32 functional. O-poor conditions are assumed. 



respectively. Our values agree well with the value of 0.8 
eV found in Ref. i and that of 0.9-1.0 eV in Ref. [13. 
Thus, our results confirm that the formation energy of 
the O vacancy in O-poor conditions is small enough to 
lead to a noticeable concentration of these defects. 

At variance with these results, Janotti and Van de 
Walle reported much higher formation energies for the 
neutral Vq.^'' They used an extrapolation procedure 
based on hDA+Ud and an additional assumption about 
the behavior of the formation energy of the charged va- 
cancy upon the band-gap correction. While the former 
extrapolation has been criticized due to the unphysical 
values to which the Ud parameter extrapolates to,— we 
argue here that it is the latter assumption that is incon- 
sistent with the hybrid functional calculations. Indeed, 
the extrapolation procedure adopted in Ref. [13 leads to 
charge transition levels that agree well with those ob- 
tained with hybrid functionals. 

The dependence of the formation energy on the elec- 
tron chemical potential is shown Fig. [2] for oxygen-poor 
conditions. For simplicity, the oxygen chemical poten- 
tial was set to the average value derived from the two 
definition schemes described above. The (-1-2/0) charge 
transition level occurs at £v + 2.38 eV. This result agrees 
well with other calculations based on hybrid function- 
als. Oba et al. found the (-1-2/0) charge transition level 
at £v + 2.23 eVr^ using the Heyd-Scuseria-Ernzerhof 
(HSE) hybrid functional based on screened exchang^ 
in which the fraction of non-local exchange was set to 
0.375 (HSE-37.5). Using the same functional but with a 
set to 0.40 (HSE-40), Clark et al. obtained the transition 
level at ey -I- 2.34 eV.'^^ Thus, it appears that when a in 
either PBEh or HSE functionals is tuned to reproduce 
the experimental band gap, one consistently obtains the 
(+2/0) charge transition level at 2.23—2.38 eV from the 
VBM. The occurrence of such an agreement has recently 
been rationalized in general termsj^ Janotti and Van de 
Walle, who adopted an extrapolation method based on 
LDA+Ud, found this charge transition level at 2.17 eV, in 
a fair agreement with the hybrid functional calculations. 

As already noted in the literature^iiii^i the charge 
transition level at ey + 2.2-2.4 eV is in stark disagree- 
ment with calculations based on other methods for cor- 



recting the band gap. For example, adopting a LDA+Ud 
schemCfSi Lany and Zunger obtained the charge transi- 
tion level at ey + 1.3 eV.^ In the LDA+C/^ method, the 
Hubbard Ud term acts on the Zn 3d states and the band- 
gap problem is not fully corrected. When one tunes the 
Ud parameter so that the position of Zn 3d states are cor- 
rectly positioned with respect to the VBM, one obtains a 
band gap of 1.5 eV, considerably smaller than the experi- 
mental one. The remaining band-gap error was corrected 
by an upward shift of the CBM.^ 

In another study, Paudel and Lambrecht adopted a 
LDA+Us/d scheme, in which the Hubbard U term was 
applied to both Zn 3d and Zn 4s states. While this 
scheme brings the theoretical band gap in agreement with 
experiment, the (-1-2/0) charge transition level is found at 
ey + 0.8 eV. Some of the results obtained in Ref. [ll| have 
recently been reviewed and improved by Boonchun and 
Lambrecht."^ We here mainly elaborate on the original 
results, but the conclusions that we draw are independent 
of this choice. Using a similar method as that of Paudel 
and Lambrecht, Lany and Zunger have obtained a level 
at ev + 0.6 eV.'' 

The charge transition levels obtained with different 
methods are compared in Fig. ^a.). We note that the 
observed differences do not stem from different electron 
densities of the defect state, as the oxygen vacancy is 
characterized by a fully symmetric state of ai symme- 
try which is well described in all schemes. The origin 
of this apparent disagreement between various methods 
has lately been discussed to some extent^ However, it 
remains unclear whether the observed differences origi- 
nate from failures of some specific methods or whether 
they point to a more fundamental problem common to 
all approximate electronic structure methods. 

A clue to the understanding why different methods 
seemingly differ so much is provided by the realization 
that the band edges of bulk ZnO calculation undergo 
drastically different shifts when going from LDA/GGA 
calculations'^^ to band- gap corrected schemes. Such shifts 
between two different electronic structure calculations 
are properly defined through the alignment of the av- 
erage electrostatic potential. For example, the LDA+Ud 
method of Ref. § yields a shift in the VBM, Aey = -0.7 
eV. The LBA+Us/d method of Ref. lll| gives a shift of 
-1-0.1 eV, while our calculations yield —1.8 eV. 

In Fig. Elb) we show the comparison of the (2-I-/0) 
charge transition level obtained with various methods, 
when the VBMs in the LDA/GGA calculations are 
aligned. This is equivalent to aligning the electrostatic 
potential of all calculations (see Sec.lIVp. With this align- 
ment, the various methods yield charge transition levels 
differing by at most 0.4 eV. This is to be contrasted to 
the variation of up to 1.8 eV achieved when the elec- 
tronic structures are aligned via their respective VBM 
[Fig.[3ja)]. Thus, these theoretical calculations do not in 
fact differ as much as has been previously claimed. Our 
conclusion is that, when a suitably defined common ref- 
erence level is adopted, the charge transition levels are 
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FIG. 3: (Color online) Calculated positions of the (+2/0) 
charge transition level of the oxygen vacancy in ZnO through 
different band-gap correction methods: (a) the various calcu- 
lations are aligned via the VBM after the band gap correc- 
tion are applied; (b) all the calculations are aligned through 
the VBM prior to shifts of the band edges required for the 
band gap correction. The illustrated results are taken from: 
"Paudel & Lambrecht '08" - Ref. fll|, "Lany & Zunger '08" - 
Ref.ll, "Lany & Zunger '07" - Ref. HT "Janotti & Van de WaUe 
'07" - Ref. [13, "Oba et al. '08" - Ref . [13, "Clark et al. '10" - 
Ref. [21I . The theoretical method is indicated in parentheses. 

more accurately described than the bulk band edges <^ In 
Sees. fVl and IVII below, we give a detailed explanation of 
this behavior and address its general consequences for 
theoretical studies of defects. 



IV. ALIGNMENT OF BULK BAND 
STRUCTURES 

The previous discussion relied on the assumption that 
the bulk band structures of two theoretical calculations 
can be aligned with respect to each other, as done in Fig. 
[Sljb) . This alignment allows one to determine the shifts in 
the valence band Aey and in the conduction band Aec 
for a given theoretical scheme with respect to another 
one. In this section, we discuss the meaning of such an 
alignment i^i^ 

The alignment between the electronic structures of the 
same bulk material within different theoretical schemes 
could in principle be achieved through the identification 



of a common reference potential. For instance, the vac- 
uum level could serve this purpose, requiring the explicit 
consideration of the surface between the considered ma- 
terial and vacuum within both theoretical schemes. Since 
the surface dipole depends on the specific crystal surface 
which is considered, the same orientation has to be cho- 
sen for both theoretical schemes. In this way, properly 
defined bulk levels in the two schemes, such as ey and £cj 
can be positioned with respect to the vacuum level and 
thus aligned. By constructon, the alignment achieved in 
this way is not an intrinsic bulk property of the two the- 
oretical schemes. Indeed, differences between the surface 
dipoles in the two surface calculations directly affect the 
alignment. 

While such a procedure can always be carried out, 
we note that the alignment between different electronic 
structures for the same bulk material is a meaningful con- 
cept only as long as their associated electron densities 
are identical (or very close). Indeed, different electron 
densities at surfaces of the material could yield different 
surface dipoles and thus the achieved alignment would 
depend on the particular surface adopted and give rise 
to ambiguity. Moreover, different surface dipoles could 
result from different electron densities in the bulk, for in- 
stance because of different theoretical equilibrium lattice 
parameters. In such a case, the alignment with respect 
to the vacuum level would again be surface dependent. 
When comparing electronic structures of bulk materials 
as achieved within different theoretical schemes, we will 
thus additionally assume that their electron densities do 
not differ essentially. In practical calculations involving 
semilocal and hybrid density functionals, this condition 
is close to being satisfied. Indeed, surface and interface 
dipoles in a variety of cases were found to differ by at 
most a few tenths of an eVi^i^i^"— 

Under the assumption of yielding close electron densi- 
ties, two different theoretical schemes can be expected to 
give similar surface dipoles. This implies that an align- 
ment to the vacuum level is equivalent to an alignment 
to the average electrostatic potentials within the bulk of 
the materials 4S This consequence is particularly conve- 
nient and allows us to compare different bulk calculations 
without the necessity of performing surface calculationsj^ 
Note, however, that it is implicitly understood that align- 
ment shifts resulting from slight differences in the elec- 
tron density are negligible when compared to the shifts 
undergone by the band edges. 

To produce Fig.[3l^b), we relied on shifts Aey and Aec 
calculated in the respective papers. Indeed, the position 
of the VBM and the CBM in the more advanced theory 
were generally given with respect to the (semi-)local den- 
sity functional calculation (LDA or GGA) for an align- 
ment with respect to the average electrostatic potential. 
For instance, the LDA+Us/d and LDA band structures 
obtained in Ref. [ill , corresponding to the left column in 
Fig. [2Ib) , were aligned through the average electrostatic 
potential in the bulk. In Refs. U^, corresponding to the 
results in the middle column in Fig. [3Jb), the authors 
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determined the shifts of the bulk bands in the LDA+Ud 
with respect to the LDA by referring the energies to O 
2s states which do not directly couple to the d states on 
which the Hubbard correction was applied. This is again 
equivalent to the alignment to the average electrostatic 
potential in the bulk. In our own calculations, presented 
in the right column in Fig.[3ljb), we aligned the two band 
structures through the average electrostatic potential in 
the bulk. Unfortunately, the reported data did not allow 
us to establish the relative alignment for all the studies 
referred to in Fig. [Sja). However, we can assume that 
similar theories yield close Aey and Aec- For instance, 
the LDA+Us/j_ calculations of Lany and Zunger— are ex- 
pected to yield similar shifts as those found by Paudel 
and LambrechtW [Fig. El^a)]. As far as the screened hy- 
brid functionals are concerned [Fig. |31Ja)], a recent study 
has shown that these functionals yield very similar shifts 
as the unscreened functionals used in our calculations, 
as long as the fraction of nonlocal exchange is tuned 
to reproduce the experimental band gapi^ Hence, al- 
though the results in Fig. ^h) are restricted to those 
studies which explicitly give the shifts in the band edges, 
the present considerations are expected to carry a much 
broader validity and to equally hold for all other calcu- 
lations reported in Fig.[3l[a). 



V. LOCALIZED AND DELOCALIZED STATES 
IN APPROXIMATE DENSITY FUNCTIONAL 
SCHEMES 

We showed above that different theoretical models give 
quite consistent results concerning the description of the 
(+2/0) charge transition level of the O vacancy in ZnO 
provided they are aligned through the average electro- 
static potential, taken as a common reference level. To 
understand why this happens, we first discuss funda- 
mental differences between localized (atomic-like) and 
extended (bulk-like) states in approximate density func- 
tional schemes. 

For (approximate) density functionals Janak's 
theorem^ applies: 



(4) 



i.e. the derivative of the total energy with respect to the 
change of occupation number /j of the highest-occupied 
state i is equal to the single-particle eigenvalue of this 
state £i, when the latter is referred to the average local 
potential."^ 

The integral form of Janak's theorem is 
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(5) 



where E^^^ is the total energy of the system with TV elec- 
trons. In the above expressions, we suppressed the spin 



variable. While in the original derivation of Janak's theo- 
rem the functionals were implicitly assumed to be contin- 
uous, Eq. dU equally applies to functionals which possess 
a discontinuity^Si^ at integer number of electrons. In this 
case one has to distinguish between left and right deriva- 
tives and the corresponding single particle eigenvalues. 
The integral form of Janak's theorem, Eq. ([S]) , applies to 
discontinuous functionals without modifications. 

In the case of localized states, such as, e.g., in 
molecules and atoms, the single particle eigenvalue in ap- 
proximate density functional schemes depends sensitively 
on the fractional occupation. Accordingly, total energy 
differences pertaining to the change of number of elec- 
trons are given by Eq. ([SJ. In particular, the ionization 
potential (IP) of a system is given by 



IP E, 



^tot - 



£N{f)df, 



(6) 



where ejv is the highest occupied orbital of the iV-electron 
system. Similarly, the electron affinity (EA) can be ex- 
pressed as 



EA = E, 



N 



E, 



N+l 



£N+i{f)df, 



(7) 



where Eat+i is the lowest unoccupied state of the A^- 
electron system. 

It has been known for some time that total energy 
differences pertaining to the change of charge state of 
a localized state are quite accurately described in ap- 
proximate density functional schemes, both in semilo- 
cal and hybrid ones. For example, Curtiss et al. cal- 
culated IPs and EAs for a large set of molecules using 
GGA and hybrid functionals^ They calculated these 
quantities via total energy differences (ASCF method), 
yielding an average deviation with respect to experiment 
lower than 0.2-0.3 eV for both GGA (BLYP) and hybrid 
(B3LYP) functionals. This accuracy is achieved despite 
the fact that the single-particle eigenvalues of the highest- 
occupied molecular orbital (HOMO) Ehomo and of the 
lowest-unoccupied molecular orbital (LUMO) elumo are 
substantially different in the GGA and in hybrid func- 
tional schemes. A similar agreement with experiment 
also holds for screened hybrid functionals <^ However, 
plain LDA yields slightly larger errors, of the order of 
0.5-0.6 cV for the same quantities. 

Wc illustrate this property in the case of the pentacene 
(C22H12) molecule in Fig. |4l^ Pentacene is a convenient 
example because, unlike several smaller acenes, it pos- 
sesses a positive electron affinity. The single-particle 
HOMO and LUMO levels, calculated with the semilocal 
PBE functional (left, solid lines), do not agree well with 
the negative of the experimental IP and EA (right). In 
particular, the single-particle gap E^^ — Elumo— Ehomo 
of 1.12 eV is severely underestimated with respect to the 
experimental gap Eg — IP — EA of 5.29 eV. The use of 
the hybrid PBEO (i.e. PBEh-25) functional (left, dashed 
lines) gives some improvement, but the calculated single- 
particle HOMO-LUMO gap of 2.34 eV remains much 
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FIG. 4: (Color online) Frontier-orbital diagram of the pen- 
tacene molecule. Left panel: HOMO and LUMO single- 
particle eigenvalues as obtained with the semilocal PBE func- 
tional (dashed, blue) and with the hybrid PBEO functional 
(solid, red). Middle panel: Ionization potentials and electron 
affinities calculated with the same functionals. Right panel: 
Experimental values for the electron affinity and the ioniza- 
tion potential. 



smaller than the experimental one. At variance, when 
calculated via total energy differences, the IPs and EAs 
in both PBE and PBEO are much closer to their corre- 
sponding experimental values, 6.64 eV~ and 1.35 eV,-- 
respectively. The two theoretical values (Fig. S]) differ 
by less than 0.10 eV, with the hybrid functional calcula- 
tion in slightly better agreement with the experimental 
results. The residual differences between calculated and 
measured values (~0.45 eV for the IP and ~0.25 eV for 
the EA) can be accounted for by the quite large elec- 
tron correlation effects in the pentacene moleculei^ In 
any case, the present result shows that these theoretical 
schemes yield total energy differences in good agreement 
with experiment and with each other, while the single- 
particle levels in the two schemes are very different. This 
is consistent with the general trend found by Curtiss et 
al. for a large set of smaller molecules 42^ 

Thus, we conclude that total energy differences per- 
taining to the change of charge state of localized states 
are accurately described with approximate density func- 
tionals. Approximating the integrals appearing in Eqs. 
© and ([7]) through the trapezoidal rule, we arrive at the 
following expressions for the IP and the EA: 



IP = E, 



N 



and 



EA = E, 



N 



E 



E, 



N+1 



-eN (I) 



-£n+i (I) 



(8) 



(9) 



Here, Sn — Ehomo and En+i — elumo- Electronic states 
at half-filling correspond to Slater-transition states.''^ 
Since Eqs. ^ and © apply equally well to various 
semilocal and hybrid functionals, the generally good 
agreement with experiment implies that the respective 
eigenvalues e(f) defined as a function of filling all ap- 
proximately cross at half-filling. This has indeed already 
been observedi^ 



The reason for this good performance of approximate 
density functionals should be ascribed to the fact that 
such functionals fulfill several exact constraints of the 
many-body fermionic system;^ In particular, the most 
relevant in this context is the generalized sum-rule of 
the exchange-correlation hole. This rule holds for sys- 
tems with an integer number of electrons, i.e. for closed 
systems in which no exchange of electrons with the en- 
vironment occursi^ This condition is enforced for most 
approximate functionals, including the LDA and various 
GGAs. Furthermore, since this constraint is naturally 
fulfilled in the Hartree-Fock theory, it also holds for any 
hybrid functional with an exchange energy of the type 
given in Eq. 

The situation is very different in the case of infinitely 
extended bulk-like states. Indeed the band-gap problem 
pertaining to (generalized) Kohn-Sham eigenvalues can- 
not be overcome by considering total-energy differences. 
When a fraction / of an electron or even a full electron 
is added to or removed from an extended state, the total 
electron density changes negligibly. Thus, the local po- 
tential, both the Hartree and the approximate exchange- 
correlation potential remain unaffected. As a result, the 
single-particle eigenvalues do not depend on the filling / 
of this state. Using the integral form of Janak's theorem 
given in Eq. ([S]), we get for the valence band maximum: 



ev - E^ot 



E, 



N-l 



and for the conduction band minimum: 



E, 



N 



(10) 



(11) 



To illustrate this property, we show in Fig. [5] the VBM 
and CBM of a-quartz calculated via total energy differ- 
ences as a function of the supercell size. The consid- 
ered cells contain 72, 144, 288, and 576 atoms, and their 
Brillouin zones are sampled at the sole P point. The 
semilocal PBE functional was used. In the case of a- 
quartz, total energy differences are very close to single 
particle eigenvalues already for the smallest cells. For 
the 72-atom cell, the difference is 0.015 eV for the VBM 
and 0.035 eV for the CBM, while for the 576 cell these 
are 0.003 eV and 0.004 eV, respectively. The particular 
case of hybrid functionals has been addressed in detail in 
Ref. [13. Hence, unlike for the localized states in Fig. SI 
the consideration of total-energy differences in the case 
of extended states is not useful to improve the compar- 
ison with experiment and the same limitations pertain- 
ing to the single-particle eigenvalues (band-gap problem) 
are encountered. ^''■^^ A similar comparison involving ex- 
tended states of GaAs and localized states of the F atom 
can be found in Ref. 4. 

The above discussion highlights an important differ- 
ence between localized and extended states as described 
within approximate density functional schemes. While 
the band-gap problem associated to single particle eigen- 
values can be circumvented by considering total-energy 
differences for localized states, such a solution does not 
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FIG. 5: (Color online) Band edges of crystalline Si02 {a- 
quartz) calculated via total energy differences as a function of 
where A^at is the total number of atoms in the supercell. 
Calculations have been performed with the semilocal PBE 
functional. 



first term can to a very good approximation be related 
to the total energy difference pertaining to a localized 
state. Formally, the first term describes the total energy 
of the vifhole manifold of states involving both defect and 
bulk states, but can be related to the localized defect 
state through the Slater transition-state approximation.— 
For atomically localized defect states, this is a very good 
approximation.^® In view of the following discussion, it 
is convenient to rewrite Eq. (1121) as 



e(+/0) = {El 



(ev - 4>) 



e(+/0) -ev. 



(13) 



where the charge transition level £(+/0) and the VBM 
Ev are referred to the average electrostatic potential (/) of 
the unperturbed bulk material. 



A. "Band-gap" problem of defect energy levels 



apply to extended states for which the band-gap problem 
remains a fundamental obstacle. Recently, a clear expla- 
nation has been put forward for justifying this different 
behavior i^iS The inaccurate total energies for large sys- 
tems with integer number of electrons stems from the 
failure of approximate density functionals in describing 
small systems with fractional charges j ^^'^'^ Indeed, ap- 
proximate functionals generally do not reproduce the 
property of the exact density functional by which the to- 
tal energy depends linearly on the number of electrons. 
There is at present an on-going effort to achieve improved 
descriptions on the basis of these ideas f^i^— The degree 
of localization required for achieving an accurate descrip- 
tion with current density functionals is still to a large 
extent an open question. We refer the reader to the in- 
teresting debate on this issue in Refs. [t^- 



VI. "THE BAND-EDGE PROBLEM" 

Having stressed the different properties of localized and 
extended states with respect to a change in electron oc- 
cupation, we return in this section to the discussion of 
charge transition levels. For the sake of simplicity, let us 
consider the (-1- /O) transition of a point defect character- 
ized by an atomically localized wave function. Using Eq. 
pO)) . we write the charge transition level e(-|-/0) as the 
difference between two terms, each of them correspond- 
ing to a total-energy difference: 



£(+/0) 



(^tot 



E, 



tot 

e: 



tot) 



( 

I ^tot, bu: 



Ik -^tot,bulk 



localized state 



dclocalized state 



(12) 



The second term clearly describes the total energy dif- 
ference pertaining to a dclocalized bulk state, while the 



Let us assume that we study the (-f/0) charge transi- 
tion level of the same defect using two different theories: 
theory I and theory II. The first theory severely underes- 
timates the band gap, while the second one gives a band 
gap in a much closer agreement with experiment. The 
two theories differ only by the exchange-correlation po- 
tential. According to Eq. (fT3|) . the corresponding charge 
transition levels referred to the respective valence band 
maxima are: 



and 



e"(+/0) = £"(+/0) 



(14) 



(15) 



We further assume that the two theories produce a suf- 
ficiently accurate representation of the electron density 
so that it is justified to align the two bulk band struc- 
tures through the average electrostatic potential <f) in the 
two theories, as discussed in Sec. HV] Under the assump- 
tion that the defect wave function differs very little in 
the two theories, we can express the difference between 
the two charge transition levels e(+/0) making use of the 
Slater transition state:'^ 



f"(+/0) - e-i(+/0) « (^d Vll - VL ^d 



(16) 



where the exchange-correlation potentials are evaluated 
with the defect state at half occupation. Only the dif- 
ference in the exchange-correlation potentials enters the 
expression in Eq. (|T6)) . Indeed, if the electron density 
and the single-particle wave functions are very similar in 
the two calculations, the interaction between the defect 
and the ionic cores, the long-range electrostatic electron- 
electron interaction, and the kinetic energy are the same 
in the two theories and cancel. 

To understand the behavior of defect levels, it is con- 
venient to focus first on defects with extremely localized 
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FIG. 6: (Color online) Schematic illustration of energy lev- 
els of various type of defect states differing by the extent of 
their wave function: (a) defect level with an atomically lo- 
calized wave function, (b) an intermediate case, and (c) an 
effective-mass-like defect. The results of two electronic struc- 
ture theories (theory I and theory II) giving different band 
gaps are compared to illustrate the band-gap problem. The 
alignment is made through the average electrostatic potential. 



wave functions. Hence, according to Eq. (|T6)) . the differ- 
ence e^^{+/0) — e^^(-|-/0) can then be expressed in terms 
of an expectation value involving the sole localized defect 
statcj^ However, we know from Sec. |V] that total energy 
differences pertaining to localized states, or, equivalently. 
Slater transition-state eigenvalues of localized states, are 
almost the same, independent of the functional. Thus, 
we get: 



(17) 



This means that charge transition levels for such defects 
are almost equal in the two theories, when the energy 
scales are aligned through the average electrostatic po- 
tential (f>. At variance, the charge transition levels are 
substantially different when the energy scales in the two 
theories are aligned through the respective valence band 
maxima, i.e. through ey in Eqs. (HH) and because 
of the different positions of the bulk band edges with re- 
spect to the potential (j). This scenario pertaining to a 
defect with an extremely localized wave function is illus- 
trated in Fig. ini by the defect state a. The validity of 
the ideal alignment illustrated by this type of defect has 
been demonstrated for a wide class of defects encompass- 
ing various host materialsi ^'^^i^^''''^' '— 

Figure M also illustrates the shifts of other type of de- 
fects. In the opposite limit, defect c corresponds to an 
effective-mass-like defect with a spatially extended wave 
function. In this case, the defect level is anchored to the 
bulk band to which it pertains and rigidly follows the 
band edge upon the opening of the band gap in theory 
II. Defect b has an intermediate extension compared to 
defects a and c, and is partially affected by the shift of 
the band edges. The relation between the departure from 
ideal alignment and the spatial extension of the defect 
wave functions has been documented for various defects 
and host materials in Ref. [^. However, the detailed be- 
havior of such defects is intrinsically system-dependent, 
and no universal considerations can be made. 

In this section, we limited the discussion to defects 
states occurring in the band gap for both theories. More 
complex situations occur when defect states are resonant 
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FIG. 7: (Color online) Schematic illustration of energy lev- 
els of various type of defect states differing by the extent of 
their wave function: (a) a deep defect level with an atomi- 
cally localized wave function, (b) an intermediate case, and 
(c) an effective-mass-like defect. The results of two electronic 
structure theories (theory I and theory II) giving the same 
band gap but different band edge positions are compared to 
illustrate the "band-edge problem". The alignment is made 
through the average electrostatic potential. 



with the band states for one of the theoriesi^ However, 
the physical description of the defect state is altered in 
such cases. The main motivation of the present work is 
to understand the effect of the band gap opening under 
the assumption that the defect wave function remains 
essentially unmodified. 



B. "Band-edge" problem for defect energy levels 

In the previous section, we compared defect charge 
transition levels as obtained within two different theories 
giving different band gaps. We found that the energy 
levels of defects states described by atomically localized 
wave functions are already well described in theories with 
a pronounced "band-gap problem" , provided those levels 
are referred to a relevant reference level. For such de- 
fects, the problem of finding the defect level is essentially 
decoupled from that of finding the band edges. 

Let us thus consider two different theories, theory I and 
theory II, yielding this time the same band gap (taken 
to coincide with the experimental one) , but giving differ- 
ent positions of the VBM and the CBM with respect to 
the average electrostatic potential (j) of the bulk. We as- 
sume that the two theories are sufficently accurate yield- 
ing in particular close electron densities, so that the en- 
ergy scales of the two theories can be aligned through ip, 
as discussed in Sec. IIVI For instance, theory I could be 
LDA+[/ in which the remaining band-gap underestima- 
tion is corrected by a rigid shift of the conduction band, 
while theory II could be a hybrid functional scheme in 
which the fraction of Fock exchange is tuned to reproduce 
the experimental band gap. For an atomically localized 
defect, the same argument holds as in the previous sec- 
tion and the charge transition levels obtained within the 
two theories are expected to fall very close to each other, 
as illustrated in Fig.[7]for defect a. In Fig. [71 a departure 
from the ideal alignment is seen for defects b and c, cor- 
responding to defect wave functions of intermediate and 
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efFective-mass-like extensions, respectively. 

Figure [7] summarizes the principal finding of the 
present work. In a condensed form, the following state- 
ment can be formulated concerning the comparison of 
charge transition levels of atomically localized defects. 
Despite the good description of the experimental band 
gap in both theories, the defect levels differ when re- 
ferred to their respective VBM, because the band edges 
in the two theories are located differently with respect 
to the common electrostatic potential (f). This occurs 
even when the defect wave function is almost identical in 
the two theories. This alignment property deteriorates 
with the extension of the defect wave function. Thus, the 
correct description of band edges relative to the average 
electrostatic potential is a crucial prerequisite for an ac- 
curate location of charge transition levels within the band 
gap. We refer to this issue as to the "band-edge prob- 
lem" for the calculation of defect levels. In other words, 
there is not only a "band gap problem" related to the 
underestimation of the band gap but also a "band-edge 
problem" related to the position of the band edges with 
respect to the average electrostatic potential, ultimately 
corresponding to an absolute alignment with respect to 
an external vacuum level. 

As far as the determination of the (+2/0) charge 
transition level of the oxygen vacancy in ZnO is con- 
cerned, the present considerations appear confirmed [cf. 
Fig. [3Kb)]. This defect level behaves like the defect state 
b in Fig. [3 showing a shift which does not depart in 
a significant way from the case of ideal alignment (de- 
fect state a). Indeed, when referred to a common ref- 
erence level, all previous calculations yield the (-1-2/0) 
level within 0.4 eV^i^"— liii^i which corresponds to just 
one ninth of the band gap of bulk ZnO. Hence, contrary 
to previous claims, we find that all previous defect cal- 
culations agree quite well with each other. In fact, these 
calculations differ in the positions of the bulk band edges 
with respect to the average electrostatic potential. 



C. Which band edge shifts are the right ones? 

These considerations lead to the question about which 
theoretical description should be adopted for position- 
ing the band edges. This corresponds to determining 
the shift Aev of the valence band and the shift Aeq of 
the conduction band, when taking the LDA or the GGA 
as a starting point. A direct comparison between the- 
ory and experiment is in principle possible. The bulk 
band structure can for instance be referred to the vac- 
uum level through a surface calculation. The VBM and 
the CBM determined in this way could then be com- 
pared with ionization potentials and electron affinities, as 
obtained by means of photoelectron and inverse photo- 
electron spectroscopy. However, such measurements are 
often shrouded by very pronounced effects associated to 
charged native defects and impurities which influence the 
electrostatics and alter the alignment. More practically, 



the validity of a given theoretical scheme can be exam- 
ined addressing band offsets at interfacesj^ Band offsets 
are well-defined quantities and can generally be measured 
through a large set of experimental techniques. The com- 
parison between theoretical and experimental band off- 
sets then allows one to determine the overall accuracy 
with which such shifts are obtained within various theo- 
retical schemes 3^^>23-86 

In the absence of experimental data, the validity of the 
shifts Aev and Aec could also be assessed by comparing 
with electronic structure calculations of higher accuracy, 
such as those based on many-body perturbation theory 
(MBPT) in the GW approximation or beyond. In- 
deed, such calculations not only provide improved rela- 
tive positions of bulk bands, but also shifts of those bands 
with respect to theoretical schemes of lower level. How- 
ever, recent work has shown that the shifts of the band 
edges with respect to the average electrostatic potential 
are more difficult to converge than relative positions of 
bandsi^ Furthermore, such shifts are sensitive to various 
levels of approximation, such as, e.g., the use of different 
models for the plasmon pole to describe the frequency de- 
pendence of the dielectric function, the inclusion of ver- 
tex corrections F, and various levels of self-consistency 
on G, W, F, and the electron wave functions i^i^i^S To 
illustrate this point, we quote a recent workj^ in which 
the relative shift of the valence band with respect to the 
overall band gap correction, i.e. Aey/ASg, was found to 
range from —0.68 to —0.42 in the case of Si02, depending 
on the level of approximation in the GW scheme. Even 
for a material as simple as Si the value of Aey / AEg as 
predicted by different GW schemes ranges from —0.75 
to +0.17 M Thus, clearly more work is needed to clarify 
these issues. A systematic study of the effects of differ- 
ent levels of approximation in MBPT on the shifts in the 
band edges is thus vital for the study of defect levels. 



VII. CONCLUSIONS 

In this work, we carried out a theoretical analysis 
of the (-1-2/0) charge transition level of the oxygen va- 
cancy in ZnO. In recent years, this defect has grown 
into a benchmark case to assess the quality of various 
advanced electronic-structure theories. Indeed, common 
approximations to density functional theory, such as the 
LDA and the various GGAs, severely underestimate the 
band gap of bulk ZnO, and the treatment at a more ad- 
vanced level thus becomes crucial even for drawing qual- 
itative conclusions. However, different advanced theoret- 
ical methods applied hitherto yielded conflicting results 
regarding the position of the defect level in the band gap. 

We here showed that apparently conflicting theoret- 
ical results are in a much better agreement with each 
other when the charge transition levels are aligned with 
respect to the average electrostatic potential rather than 
to the respective valence band maximum. We showed 
that the former alignment is equivalent to the choice of 
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a common external potential such as the vacuum level, 
provided the electron densities are sufficiently accurately 
described. We have rationalized our finding by consid- 
ering fundamental differences between the ways approx- 
imate density functionals describe localized and delocal- 
ized states. For localized states, the "band-gap problem" 
can generally be overcome through the consideration of 
total energy differences. On the other hand, such a solu- 
tion is not applicable to dclocalizcd states, for which the 
"band-gap problem" remains an intrinsic shortcoming. 

In particular, the present study highlights a specific 
criterion that needs to be fulfilled in order to properly 
describe charge transition level and formation energies of 
defects. We clearly demonstrated that the band struc- 
ture of the host material needs to be correctly positioned 
with respect to an external potential, such as the vacuum 
level. When the electron density is accurately describcid, 
this alignment condition can equivalently be replaced by 
the alignment with respect to the average electrostatic 



potential in the bulk. This condition is additional with 
respect to the accurate reproduction of the band gap. 
Our analysis of the oxygen vacancy in ZnO shows that 
conflicting theoretical results arise for theories yielding an 
accurate band gap, but differing positions for the band 
edges. 
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